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(57) Abstract: A process for performing an improved mixed fre- 
quency-time algorithm to simulate responses of a circuit that re- 
ceives a periodic sample signal and at least one information sig- 
nal, the process selecting a set of evenly spaced distinct time points 
(804) and a set of reference time points (806). Each of the refer- 
ence points is associated with a distinct time point, and a reference 
time point is a signal period away from its respective distinct time 
point. The process finds a first set of relationships (808) between 
the values at the distinct time points and the values at the refer- 
ence time points, and a second set of relationships (810) between 
the values at the distinct time points and the values at the reference 
time points. The first and second set of relationships are combined 
to establish a system of nonlinear equations in terms of the values 
at the distinct times only (812). The nonlinear equations are solved 
to find simulated responses of the circuit in the time domain, and 
the simulated circuit responses are then converted to a frequency 
domain (814). 
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METHOD AND APPARATUS FOR SIMULATING 
QUASI-PERIODIC CIRCUIT OPERATING CONDITIONS USING A 
MIXED FREQUENCY/TIME ALGORITHM 

5 FIELD OF THE INVENTION 

The present invention relates generally to analog circuit design simulations, and more 
specifically to analog circuit design simulations using a mixed frequency/time approach. 

BACKGROUND OF THE INVENTION 

Using a description language such as a netlist and device models an analog circuit can be 

10 first designed in terms of its predetermined inputs and expected outputs. The analog circuit 
design is then simulated before it is physically fabricated on a silicon chip. 

One of the most difficult challenges in analog circuit simulation is the analysis of the 
circuits that operate on multiple time scales. Typical examples of this type of circuits are 
switched-capacitor filters and circuits used in RF (radio frequency) communications systems. 

1 5 Applying standard transient analysis to a circuit of this type requires simulation of the detailed 
responses of the circuit over hundreds of thousands of clock cycles (millions of time points). 

Many circuits of engineering interest are designed to operate near a time-varying, but 
quasi-periodic, operating point. Some of these circuits can be analyzed under the assumption 
that one of the circuit inputs produces a periodic response that can be directly calculated by 

20 steady-state algorithms, thus avoiding long transient simulation times. Under this assumption, 
all other time-varying circuit inputs are treated as small-signal by linearizing the circuit around 
the periodic operating point. 

Existing algorithms are able to find periodic operating points and to perform periodic 
time-varying small-signal analysis. However, many circuits cannot be analyzed with the 

25 periodic-operating-point-plus-small-signal approach, because the above-described assumption 
may not apply. For example, predicting intermodulation distortion of a narrowband circuit, such 
as a mixer-plus-filter circuit, involves calculating the nonlinear response of the mixer circuit, 
driven by an LO (local oscillator), to two high-frequency inputs that are closely spaced in 
frequency. The steady-state response of such a circuit is quasi-periodic. 

30 The analog circuit simulation is farther complicated by the fact that many multi-timescale 

circuits have a response (again mixers and switched-capacitor filters are typical examples) that is 
highly nonlinear with respect to at least one of the exciting inputs, and so steady-state 
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approaches, such as the multi-frequency harmonic balance approach, do not perform well. To 
circumvent these difficulties, mixed frequency-time (MFT) algorithms have been proposed. 
Specifically, the MFT algorithms exploit the fact that many circuits of engineering interest have 
a strongly nonlinear response to only one inout, such as the clock in the case of a switched- 
5 capacitor circuit, or local oscillator in the case of a mixer, but respond only in a weakly nonlinear 
manner to other inputs. 

Unfortunately, existing MFT algorithms suffer from several drawbacks that prevent their 
application to practical circuits, particularly large circuits. In existing MFT algorithms, poor 
sample point selection leads to ill-conditioned simulation environment, in which simulation 
1 0 values may be unsolveable with acceptable accuracy. In addition, existing MFT algorithms are 
based on a matrix-explicit linear solver (via Gaussian elimination) whose computational cost (or 
time) is proportional to an order of N 3 for each Newton iteration, where N is the number of nodes 
of the circuit in simulation. 

A new class of algorithms has been developed for simulating multi-timescale circuits by 
15 converting the circuit DAE (differential-algebraic equation) into an equivalent multi-variable 
partial differential equations (M-PDE). However, the effectiveness of the M-PDE method to 
simulate large circuits has yet to be proven. In addition, there is evidence that, for some circuits, 
the M-PDE method generates inaccurate simulation results. 

There is, therefore, a need in the art for a method and apparatus that utilizes the MFT 
20 method to accurately simulate large circuits. 

There is another need in the art for a method and apparatus utilizing the MFT method to 
simulate large circuits with reduced computational cost and increased speed. 

There is still another need in the art for a method and apparatus for generating an efficient 
linear problem solver structured such that the MFT method can accurately simulate large circuits 
25 with improved convergence. 

The present invention provides a method and apparatus to meet these and other needs. 

SUMMARY OF THE INVENTION 

To overcome the shortcomings of the available art, the present invention discloses a 
30 novel method and apparatus for simulating analog circuits by using a mixed frequency/time 
approach. 

In broad terms, the present invention provides a method for simulating responses of a 
circuit, the circuit receiving a periodic sample signal and at least one information signal. The 
method comprises the steps of: selecting a set of distinct time points; defining a set of reference 
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time points, wherein each of the reference time points is associated with one of the distinct time 
points; establishing a first set of relationships between the values at the distinct time points and 
the values at the reference time points; establishing a second set of relationships between the 
values at the distinct time points and the values at the reference points; combining the first and 
5 second relationships to establish a system of equations in terms of the values at the distinct time 
points; and finding responses of the circuit at the distinct time points by solving the established 
system of equations. 

The present invention also provides a corresponding apparatus for performing steps in the 
method described above. 

10 

BRIEF DESCRIPTION OF THE DRAWING 

The above mentioned advantages of the present invention as well as additional 
advantages will be more clearly understood as a result of a detailed description of the preferred 
embodiments of the invention when taken in conjunction with the following drawing in which: 
15 Figure 1 shows circuit 102, which will be used to illustrate a simulation process in 

accordance with the present invention; 

Figure 2 shows an exemplary envelope in response to two inputs shown in Figure 1 ; 

Figure 3 shows a scheme of selecting K (K = 2) sample points along the sample envelope 
shown in Figure 2 to obtain five distinct time points, in accordance with the present invention; 
20 Figure 4 shows the effectiveness of the preconditioning process in compressing the 

eigenvalues for an RF receiver circuit, in accordance with the present invention; 

Figure 5 shows the effectiveness of the preconditioning process in reducing the number 
of iterations needed to solve each MFT Newton update equation, in accordance with the present 
invention; 

25 Figure 6 shows a simulation result for a filter circuit in accordance with the present 

invention; 

Figure 7 shows the intermodulation distortion of a high-performance receiver; 

Figure 8 shows a flowchart illustrating a process of simulation the responses of a circuit, 
in accordance with the present invention; 
30 Figure 9 shows a flowchart illustrating an exemplary process of solving the system of 

equations established in Figure 9, in accordance with the present invention; and 

Figure 10 shows block diagram illustrating an exemplary computer system, which can be 
used as a hardware platform for executing the program that performs the processes shown in 
Figures 8 and 9, in accordance with the present invention. 

3 
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DETAILED DESC RIPTION OF PREFERRED EMBODIMENTS 

Figure 1 shows circuit block 102 including N circuit nodes, which will be used to 
illustrate a simulation process in accordance with the present invention. Circuit 102 receives S 
information signals II, 12, Is from inputs 104.1, 104.2, 104.s, respectively. Circuit 102 
5 also receives a clock (or reference) signal C from input 106. In response to receipt of the input 
signals, circuit 102 generates a signal at output 108. 

1. The improved MFT algorithm 

The behavior of a circuit (such as circuit 102 in Figure 1) can be described by a set of 
1 0 nonlinear differential-algebraic equations (DAEs) that can be written as 

^fi(v(0) + /(v(0) + fi(0 = 0, (1) 

Where Q(v(t)) e 9t* is typically the vector of sums of capacitor charges at each node, 

I(v(t)) e 9t N is the vector of sums of resistive currents at each node, u(t) e SR* is the vector of 

1 5 inputs, v(t) e 9* N is the vector of node voltages, and N is the number of circuit nodes. 

The present invention is particularly advantageous in situations where the input signal 
lift) is quasiperiodic. A signal is Z-quasiperiodic if it can be written as a Fourier series with L 
fundamental frequencies. RF circuits are generally influenced by one periodic timing signal, 
often referred to as the LO (local oscillator) or the clock, and one or more information signals. If 

20 / c denotes the clock signal frequency (such as the signal at input 106 shown in Figure 1), and 

fx . . .f s are the S information signals (shown in Figure 1), then the (S + l>quasiperiodic input can 
be written as: 

25 • * — m 

In a preferred embodiment, the present invention utilizes two conditions to improve the 
simulation of quasi-periodic circuit operating conditions. The first condition is that the circuit of 
30 interest possesses a quasiperiodic steady-state response. That is, v(t) is an S + 1 quasiperiodic 
signal with fundamentals f\ ,. . .f s ,f c .. The second condition is that all physical circuits have a 
finite bandwidth. Using these two conditions, the present invention selects only a finite number 
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of Fourier series terms to approximate v(t) while maintaining the necessary accuracy. Thus: 

K K 

v(')= £ - E Z V{k x ,...,k„k e ) (3) 

x e -/2^,/,l... e -J2zk e f e t 

where fe.AJ e C*. (An interesting property of the MFT algorithm is that it is not 

necessary to truncate to a finite number of harmonics of f c ) 

Assume that v(t) is sampled at a discrete set of points /'„ = t 0 + nT Ct where T c = l// c is the 
clock period, t 0 e [0,r c ) and n runs over the integers, to obtain a discrete signal v(t% Since 



K K 

r ('' n )= X - 2 v(K-,k s )x e-J 2 **'* 1 - ■■■e- J2 " kS/s " < 4 > 



where 



1 5 the "envelope" v(0 is S-quasiperiodic and can be represented as a Fourier series in only the 
"information" fundamentals. The clock fundamental has disappeared, thereby forming an 
envelope. 

Figure 2 shows an exemplary envelop in response to the signals at input 104.1 (assuming 
that the inputs at 104.2, 104.S are zero inputs) and input 106. In Figure 2, the hollow, 
20 circular dots represent samples (or distinct points), while the continuous waveform (dashed line) 
is the waveform having V as its Fourier coefficients, or equivalently, obtained by Fourier 
interpolation of the sample points. 

In principle, because there are only K = tf s= j (2K 3 + 1) Fourier coefficients to represent 

v, once the value of K distinct points t lt ... ,t k along the sample envelope are known, then the full 
25 envelope can be recovered. The envelope corresponding to the quasiperiodic operating point is 

obtained by obtaining K sample points that lie on the solution to the DAE given by equation (1). 
Figure 3 shows a particular scheme of selecting K(K=2) sample points along the sample 

envelope shown in Figure 3 to obtain five distinct points (or distinct time points) [5 = (Kx2) + 

1], in accordance with the present invention. In Figure 3, the distinct points are denoted by 
30 circles, and the reference points (reference time points) are denoted by squares. 
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Let us define the state transition function <p (v 0 , t h tj = v(tj) = v(tj): v(t) that satisfies 
equation (1) for t e ft h tjj and v(t0 = v 0 . In particular, define the vector 



v 0 = [v (/, ),..., v T {t K )] 7 = [v T (/, ),..., v T {t K )} T , (6) 

5 

where superscript T denotes matrix transpose, to contain v at the K sample points t k = t 0 + n h T c , 
k= l...K,n k eZ. The value of the K points that follow by one cycle can be obtained from the 
transition function, 

ltfT c =^ T ^ + 7 U<... ,v'(, A . +T C )] T =[ < *(v(/ 1 ),/ l ,r l + r c ) r ,... J(v(t K ), t K + T C ) T ] T 

(7) 



which may be written more compactly by introducing the multi-cycle transition function that is 
the collection of the K transition functions from tk to tk + T a , as 

v r e = * n (v 0 ). (8) 

1 5 Note that for each mode n, the vector of signals on that node, at the sample time plus one 

clock cycle, v T , is a delayed version of the signals at the sample points (this will be discussed in 
greater detail below). There exists, therefore, a linear operator D Tc that maps to v£ 



v£ = D T v%. 

(9) 

20 

Note that D Te is a real matrix and independent of node n. Hence equation (9) holds for 
each n = 1,. . N, and represents a boundary condition on a solution to (1). 
Combining equations (8) and (9) gives 



25 (D Te ® /„)v 0 -<J> r (v o ) = 0, (10) 



6 
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where ® is the Kronecker product 1 and 1 N is the N by N identity matrix. Equation (10) is a 
system of KN nonlinear equations and KN unknowns that can be solved for the envelope sample 
points. From these sample points and the transition functions, the circuit's quasiperiodic 
operating point (in particular, the spectrum of v) can be recovered. 

2. Sample point selection 

To construct the matrix D Te , referred to as the delay matrix, consider the Fourier series of 
vo and vt c . Referring to equation (4), equation (11) holds 

A* A' ' 

v(/'„+r c )= £ ... jT F(*,,... ,k s )e-W'... e-^'-Q^t, * s )(ll) 



A: t = — AT j k s =>— K 



where 

Thus if r is the matrix mapping sample points on the envelope to Fourier coefficients, 
then the delay matrix may be constructed as 

D Te = r _, Q r r. (13) 

In particular T may be constructed as the Kronecker product of one-dimensional (2^+7)-point 
Fourier-transform matrices 

Y ( J ) = e J 2ftm f*^^ 2K ^) (14) 

as 

20 r = r (1 >®... ®r (S) < 15 ) 

From tlie properties of Kronecker products, T 1 is likewise a Kronecker product of the inverses of 
the r®. In the existing MFT algorithms, no particular consideration was given to the choice of 
the sample points t k , so that the T^s there are ill-conditioned matrices corresponding to an 
"almost-periodic" Fourier transform. By contrast, the improved MFT algorithm of the present 
25 invention performs a process of choosing well-conditioned sample points. 



7 
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Assume the K sample points can be arranged into an S-dimensional array x(k if ...,k s \ 
K s £ k s <> K Sy 1 < s << S, such that for a given dimension s, there exists an integer p, and 



r(- +l,..-)-r(- ,*,,•••) = 



2K S + 1 



+ PT, 



(16) 



5 holds. In this case, the entries of the r* J * matrices are: 



That is, they are the DFT matrices, and the matrix T : C 2K ' + 1 x - - * x C 2K ' +1 
=»C 2K ' + 1 x — x C 2K « +1 represents an S-dimensional DFT. Thus T has a condition number of 
10 one; it is perfectly well-conditioned. 

3. Matrix-implicit solution procedure 

The Newton's method can now be employed to solve 



U>r, ® /*)v 0 -*r,(v.) = 0, 



(18) 



15 At iteration z, the Jacobian matrix is given by 



d<£> 

D rr ®I N v 0 -^L 



v 0 =?4 



(19) 



Recall from (13) D Tc = r'Q^ which is fixed through all Newton iterations. Let 
J = 1 v 0 = v' 0 be obtained from the multicycle transition function by 



20 



dv 0 



v 0 ('i )=v 0 (/,)' 



dv 0 (t K ) 



7o(fx)=7o(f A )' 



(20) 
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Note that J is block-diagonal. Defining b = -(Dj c ® In)v 1 0 - <Kv$, the Newton iteration is 
performed by solving the equation 



((D Tr <8>J„)- J)A(vi) = i, 



(21) 



10 



using an iterative Generalized Minimal Residual (GMRES) solver, and setting 



— t+i _ — i 



(22) 



Each iteration of GMRES requires a matrix-vector multiplication. For a vector q e 9f at t the 
term (D Te ® I N )q is calculated by first applying a K dimensional DFT N times, then scaling each 
row with Q Te , and finally applying a K dimensional inverse DFT N times. 

Let q be partitioned into q = [q] y . . .,^] T ,qjt e <K N , for 1 < < K. Then 



dv 0 (t K ) 



qK 



(23) 



The calculation of each ^-^ can be carried through matrix-vector multiplication and 
backsolving without explicitly forming the ]f^Ji k matrix. 

4. Preconditioning 

1 5 For many problems, the GMRES algorithm is not efficient for solving equation (21) 

without an effective preconditioned To analyze the reason, consider the case where the state 
transition function of the circuit, over one clock cycle, is approximately linear, that is 
<KX>^ + Tc) ~ H%(t). Linear circuits are an obvious example of a case where this is true, and 
while nonlinear circuits will have nonlinear state-transition functions, if the method performs 

20 poorly for linear circuits it surely will not work well for nonlinear circuits either. However, 
many nonlinear circuits have a state-transition function that is nearly linear, a fact which is 
exploited below to construct an effective preconditioner. The convergence of the GMRES 
method will depend on the location of the eigenvalues of the Jacobian matrix, D T -J. If Xh is an 



BNSDOCID: <WO 01B8830A1_I_> 



WO 01/8883(1 



PCT/US01/I5185 



10 



15 



eigenvalue of the matrix H, then j°> xc -X H , where co = 2ji(/t/, + /fc/ 2 + ■ "k/ t ) will be an 
eigenvalue of D Tc - J, for every ki,k 2 , — in the MFT analysis. Thus unless all the secondary 
input frequencies are nearly commensurate with the clock frequency, the eigenvalues of D Tt - J 
will be "fanned out" by delay matrix. This will cause severe convergence problems for the 
GMRES solver. Roughly speaking, the GMRES algorithm in the MFT algorithm with A: total 
harmonics will take K times as many iterations to coverage than the GMRES iteration for the 
steady-state problem with only the clock excitation applied. This follows because the 
eigenvalues of H are typically inside the unit circle of the complex plane. The delay matrix 
replicates the eigenvalue structure # times, each shift being a complex number of order unity, 
and generally causing the convex hull of the eigenvalues of D Tc -J to enclose the origin. 

The following lemmas about the properties of Kronecker products are needed to perform 
the formal analysis. 

Lemma 5.1 If Al , M Ap e p"*" , Bl , B2 Bp e F m then 

A 1 A 2 ...A p )®(B 1 B 2 ...B p ) = (A l ®B l XA 2 ®B 2 )...(A p ®B p ). 

Lemma 5.2 IfAe F ' m , B e /"then 

(a) {A ® I n )(/„ ®B)~ (/„, ® B\A ® /„ ). 

(b) (A®By x =XA~ 1 ®B- i ). 

Theorem 5.3 If X H is an eigenvalue of dQ>ldv g , then e" 2 ™* 1 * is an eigenvalue of the MFT 
Jacobian matrix. 

The proof is as follows. For linear circuits, the diagonal blocks of J=- are the same, i.e., 
a v 0 (r,)' Denote a diagonal block as H, then the Jacobian matrix is equal to 



(r-'n^r)®^-^®^) 
= (r-' ®/ A ,)(Q 7 . ®J fl )(r®i N )-(i K ® H ) 
= (r-'®i„)w rc ®r„)-(r- t ®i N )-\i K ®H)(r®i N r i Kr®r N ) 
= (r-' ®/„){(Q n ®i„)-(r®i N )(i K ®H)(r®i N y 1 }(r®i N ) 
=(r-' ®i N )m Tt ®i„)-(i K ®Hxr®i N )(r®r N y l }(r®i N ) 

= (T-' ®/ w ){(Q 7 . ®I„)-{I K ®H)}(T®I N ). 



(24) 
(25) 
(26) 

(27) 
(28) 
<29) 



10 



BNSOOCID: <WO 018B830A1 J_> 



WO 01/88830 



PCT/US01/15185 



Equation (24) to equation (25) holds because of In = In In In and Lemma 1 . Equation (26) to 
equation (27) holds due to Lemma 2(b), and equation (27) to equation (28) holds due to Lemma 
5.2(a). Since (T l ® In) is unitary and its inverse is (r ® In)'\ the right hand side of equation 
(29) has the same spectrum as (Qj; <8>In)-(Ik® H)- It is easy to verify that 
5 (£l Te ®In)-( Ik ®H)is block diagonal, hence its eigenvalues are the union of eigenvalues of all 
the blocks, J 2 ""** I N - H, for k = 7, • • • ,K. 

The preceding analysis suggests a good way of preconditioning for solving the Newton 
equation (21). Solving equation (21) is equivalent to solving 

{ci Tt - «r ® /„) mxt 1 ® i N ))}y = n>, (3 o) 

10 where y = TA v K A good choice of preconditioner is P = (Q Te ® ./#) - ( Ik ® H) where H can be 
chosen as the Jacobian matrix from the steady-state analysis in the initial guess stage discussed 
in Section 5, or any of the diagonal blocks, ^J^, for i = 1, . . K 7 of . In particular, if the 

single-cycle state-transition function is linear and time invariant, then the Newton equation can 
be solved in a single GMRES iteration. Note that the preconditioner presented here is effective 

15 if the Jacobian of the state-transition function is nearly constant over multiple cycles. The circuit 
behavior inside each clock cycle is hidden from the preconditioner. This is not the case in, for 
example, the time- or frequency-averaged preconditioners typically used in modem harmonic 
balance codes. For this reason the preconditioner presented here may perform well under much 
weaker assumptions about the circuit behavior, in particular at higher power levels. 

20 For each GMRES iteration, a system Pu - v has to be solved. Since P is block diagonal, 

it needs to solve a sequence of K systems (e! 2 ™ k7c I N - H)Uk=Vk> for £ = 1 , . . . ,K. where 
u T =[u T /9 . . ,,w£] r and v T -[v T j9 . . -,v£] r . The preconditioner can be applied very efficiently by 
incorporating a Krylov subspace reuse algorithm, as the linear systems to be solved are the same 
as arise in the small-signal analysis for periodic time-varying systems. The basic idea of the 

25 algorithm is that the Krylov subspaces associated with the matrices e mTc - H are very similar for 
different ©k. Essentially, the Krylov-subspace re-use algorithm allows the preconditioner for the 
matrix (D Tc ® I N ) -J to be applied with only slightly more cost than an iterative solve with the 
matrix H. 

Figure 4 shows the effectiveness of the preconditioning process in compressing the 
30 eigenvalues for an RF receiver circuit. Specifically, Figure 4 shows eigenvalue distribution 
before and after the preconditioning process. In figure 4, the eigenvalues are very tightly 

11 
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clustered around unity, indicating excellent performance of the preconditioner and very rapid 
GMRES convergence. 

Figure 5 shows the effectiveness of the preconditioning process in reducing the number 
of GMRES iterations needed to solve each MT Newton update equation, for the same RF circuit 
mentioned above. Only three iterations are needed to reduce the residual by a factor of 10* 2 , 
whereas without the preconditioning process, over 400 iterations are needed to achieve any 
reduction in the residual at all. Since the MFT circuit equations are not solved exactly, on 
average there is a performance advantage in the MFT method to using approximate solves of the 
Newton update equation, and therefore GMRES is converged to a relatively loose tolerance. 

5. Improving Newton convergence 

Rapid convergence of Newton's method can only be assured with a good initial 
estimation. To achieve a good initial estimation, the present invention first calculates the 
periodic steady state response of the circuit with the clock signal applied, while suppressing 
other non-DC signals. Using the steady state solution as an operating point, a small-signal 
analysis is performed by treating non-clock fundamentals as small signals. As a result of the 
small signal analysis, amplitudes atf s + kj ct for -K S <>K S ^^,l^s^S,are generated. These 
amplitudes are transformed into time domain initial conditions via inverse multidimensional 
discrete Fourier transform (DFT). At higher input power levels, using a Newton continuation 
method, with the amplitude of the non-clock signals as the continuation parameter, is generally 
effective in securing convergence. 

6. Spectrum calculation 

After the solution is converged, the values v = [v (t if , vfc) r , " " ' 9 v(tx) T ] T and the 
integration solution in v = [v (t\f , vifif ' 'Mfj^f are available. From these pieces of 
information, the spectrum v(t) can be obtained. Let 




k s =-Ks 



(31) 



Define v(r) = [v(t { +r) r , v(t 2 +r) T 



.,v(t K +r) T ] r . Then 



v( T ) = {(r-'Q(r)) ® I„ } x Z£= -K t V{k x k s ,k c 



(32) 
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Then for each AT TV- vector V(\kc) where -K c <k c £K c which is collection of all N-vectors 
V(k h . . .,k S) kc) t where -Kj £k]Z K it ,-K s <.k s <,K, (the actual order is determined by the Fourier 
transform), 

5 r(-.* t ) = l f{(fi(r)- , n®/ w }v(V 2 "** / " (33) 

Forming {(fi (t^T)®/// } v (x) requires the values for vf//+r),. . . ,v(frfr) , or synchronized time 
steps between cycles. The total computational costs is one AW-vector integration and M Fourier 
transforms, where Mis the number of synchronized time points. 
1 0 The synchronized time step requirement may not be easily met in practice. One 

alternative is to use interpolation schemes. However, these schemes potentially lose accuracy. 
Another alternative is to utilize integration instead of multidimensional discrete Fourier 
transforms. Specifically, it is easy to verify that 



(34) 



where E p is a KNxN block matrix whose pfh NxN block is In and other blocks zero, andjp is 
determined by (k\ 9 . . .,£$) from the Fourier transform. Calculating equation (34) does not require 
synchronized time points. The total cost of calculating V(.,k) is K 2£2V-vector integrations plus 
20 one final Fourier transform. However, it might be more expensive since integrations normally 
cost more than Fourier transforms. 

7. Simulation Utilizing A Preferred Embodiment Of MFT Method 

The first example is a low-pass switched-capacitor filter of 4kHz bandwidth and having 
25 238 nodes, resulting in 337 equations. To analyze this circuit, the improved MFT of the present 
invention analysis was performed with an 8-phase 100kHz clock and a 1 V sinusoidal input at 
100Hz. 
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The 1000 to 1 clock to signal ratio makes this circuit difficult for traditional circuit 
simulators to analyze. In the improved MFT method, three harmonics were used to model the 
input signal. The eight-phase clock resulted in the need to use about 1250 timepoints in each 
transient integration. This brings the total number of variables solved by the analysis to slightly 
5 less than three million (337 x(2x3 + l)x 1250 = 2,948,750). The simulation took a little less 
than 20 minutes CPU time to finish, on a Sun UltraSparcl workstation with 128 Megabyte 
memory and a 1 67MHz CPU clock. Figure 6 shows the output spectrum of the filter. 

The second example is a high-performance image rejection receiver. It consists of a low- 
noise amplifier, a splitting network, two double-balanced mixers, and two broadband Hilbert 

10 transform output filters combined with a summing network that is used to suppress the undesired 
side-band. A limiter in the LO path is used for controlling the amplitude of the LO. It is a rather 
large RF circuit that contains 167 bipolar transistors and uses 378 nodes. TTiis circuit generates 
987 equations in the simulator. 

To determine the intcrmodulation distortion characteristics, the circuit was driven by a 

15 780MHz LO and two 50mV closely placed RF inputs, at 840MHz and 840MHz+10KHz, 

respectively. Three harmonics wee used to model each of the RF signals. 200 time points were 
used in each transient clock-cycle integration, considered to be conservative in terms of accuracy 
for this circuit. As a result, nearly ten million unknowns (987 x (2 x 3 + l) 2 x 200 = 9,672,600) 
were generated. It took 55 CPU minutes to finish on a Sun UltraSparelO workstation with 128 

20 Megabytes of physical memory and a 300MHz CPU clock. Figure 7 shows 3rd and 5th order 
distortion products. 

To understand the efficiency of the improved MFT method of the present invention, 
consider that traditional transient analysis would need at least 80,000 cycles of the LO to 
compute the distortion, a simulation time of over two days. In contrast, the MFT method of the 
25 present invention is able to resolve very small signal levels, such as the 5th order distortion 
products show in Figure 7. 

Solving the MFT equations by direct factorization methods is also impractical, as the 
storage needed for the factored rank-50,000 (987 x(2x3rl) 2 = 48,363) MFT Jacobian of 
Equation 19 is several gigabytes. Forming the Jacobian matrix by direct methods would also 
30 require computation time proportional to the cost of 50,000 transient integration cycles, again a 
number on the order of days. 

Figure 8 shows a flowchart illustrating a process of simulating the responses of a circuit 
such as the circuit shown in Figure 1, in accordance with the a preferred embodiment of the 
present invention. 
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Step 804 selects a set of evenly spaced distinct time points shown as the circle dots in 
Figure 3. The details of step 804 can be found inequation (16) and related descriptions. 

Step 806 defines a set of reference time points shown as the square dots in Figure 3. As 
shown in Figure 3, each of the reference time points is associated with a respective distinct time 
5 point, and each of the reference time points is one signal period (or clock cycle) away from its 
respective distinct time point. 

Step 808 establishes a first set of relationships between the values at the distinct time 
points and the values at the reference time points. The details of step 808 can be found in 
equation (8) and related descriptions. 

10 Step 810 establishes a second set of relationships between the values at the distinct time 

points and the values at the reference time points. The details of step 810 can be found in 
equation (9) and related descriptions. 

Step 812 combines the first and second sets of relationships to establish a system and 
equations that contain the values at the distinct time points only. The details of step 808 can be 
1 5 found in equations ( 1 0) and ( 1 8) and related descriptions. 

Step 814 finds (or generates) the simulated responses of the circuit at the distinct time 
points by solving the established system of equations. If a circuit includes N internal circuit 
nodes and M outputs, step 814 can find (or generate) the simulated responses for all of the N 
internal circuit nodes and M outputs. The details of step 814 can be found in equations (18>(22) 
20 and (30) and related descriptions. 

Figure 9 shows a flowchart illustrating an exemplary process of solving the established 
system of equations in step 814 of Figure 8, in accordance with the present invention. 

Step 904 selects a set of estimated values to reflect estimated circuit responses at the 
distinct time points. The details of step 904 can be found in Section 5. 

25 Step 906 establishes a system of linear equations at the estimated values. The details of 

step 906 can be found in equation (21) and related descriptions. 

Step 908 preconditions the system of linear equations to improve the convergence of 
solution to the system of linear equations. The details of step 908 can be found Section 4. 

Step 910 solves the system of linear equations to generate the correction values to adjust 
30 the estimated circuit responses at the distinct time points. The details of step 910 can be found in 
equations (21)-{22) and related descriptions. 
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Step 912 adjusts the estimated values as newly estimated values to reflect the estimated 
circuit responses at the distinct time points. The details of step 912 can be found in equations 
(21>-(22) and related descriptions. 

Step 914 determines whether the adjusted estimated values have an acceptable accuracy 
5 to represent the circuit responses. If the determination is negative, the process is led to step 906. 
If the determination is positive, the process is led to step 916. The estimated values and adjusted 
estimated values are in time domain. 

Step 916 converts the estimated values from time domain to frequency domain. The 
details of step 916 can be found in Section 6, equations (31) - (34). 

10 8. Hardware platform 

Figure 10 shows a block diagram illustrating an exemplary computer system 1000 ? which 
can be used as a hardware platform for executing the program that performs the processes shown 
in Figures 8 and 9. 

As shown in Figure 10, computer system 1000 includes system bus 1001, processing unit 
15 1002, memory device 1004, disk drive interface 1006, hard disk 1008, display interface 1010, 
display monitor 1012, bus interface 1014, mouse 1016, keyboard 1018, and network 
communication interface 1020. 

The hard disk 1008 is coupled to disk drive interface 1006; monitor display 1012 is 
coupled to display interface 1010; and mouse 1016 and keyboard 1018 are coupled to bus 
20 interface 1014. Coupled to system bus 1 00 1 are processing unit 1 002, memory device 1 004, disk 
drive interface 1006, display interface 1010, and network communication interface 1020. 

The memory device 1004 stores data and programs. Operating together with disk drive 
interface 1006, hard disk 1008 also stores data and programs. However, memory device 1004 
has faster access speed than hard disk 1008, while hard disk 1008 has higher capacity than 
25 memory device 1004. 

Operating together with the display interface 1010, display monitor 1012 provides visual 
interfaces between the programs being executed and users, and displays the outputs generated by 
the programs. 

Operating together with bus interface 1014, mouse 1016 and keyboard 1018 provide 
30 inputs to computer system 1000. 

The network communication interface 1020 provides an interface between computer 
system 1000 and network 104 in accordance with predetermined networking protocols. 
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The processing unit 1002, which may include more than one processor, controls the 
operations of computer system 1000 by executing the programs stored in memory device 1004 
and hard disk 1008. The processing unit also controls the transmissions of data and programs 
between memory device 1004 and hard disk 1008. 
5 In the present invention, the program for performing the steps shown in Figures 8 and 9 

can be stored in memory device 1014 or hard disk 1018. The program can be executed by 
processing unit 1002. 

9. Summary 

The present invention improves the existing MFT method. The MFT method of the 

10 present invention is an efficient approach to analyzing multi-frequency nonlinear effects such as 
intermodulation distortion. Making the MFT method computationally efficient on problems of 
engineering interest required careful construction of the delay matrix, matrix-implicit Krylov 
subspace iterative linear solvers, and a preconditioner tailored to the MFT method and the 
circuits it typically analyzes. As a result, nonlinear systems comprising tens of millions of 

15 unknowns can be solved in less than an hour with computational resources commonly available 
to engineering designers. 

One salient advantage of the MFT method in the present invention is in computing the 
functions 0> and the product of the Jacobian of O with some vectors. Both computations are 
essentially the solution of an initial value problem. Each application of the operator D Tp - Jot 

20 calculation of the Newton residual, involves solving K such initial value problems, that is, 
integrating K sets of DAEs forward in time over one clock period. Each of the £ problems, 
however, is essentially decoupled. Parallel implementations of the MFT will therefore enjoy 
very efficient processor utilization. This decoupling also assists the implementation of out-of- 
core solvers. In fact, it has been observed it is possible to implement the MFT algorithm as an 

25 out-of-core algorithm with over 80% average CPU utilization. 

While the invention has been illustrated and described in detail in the drawing and 
foregoing description, it should be understood that the invention may be implemented through 
alternative embodiments within the spirit of the present invention. Thus, the scope of the present 
invention is not intended to be limited to the illustration in this specification, but is to be defined 

30 by the appended claims. 
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1 . A method for simulating responses of a circuit, the circuit receiving a periodic 
sample signal and at least one information signal, the method comprising the steps of: 

5 (a) selecting a set of evenly spaced distinct time points; 

(b) defining a set of reference time points, wherein each of the reference time 
points is associated with one of the distinct time points; 

(c) establishing a first set of relationships between the values at the distinct time 
points and the values at the reference time points; 

10 (d) establishing a second set of relationships between the values at the distinct 

time points and the values at the reference points; 

(e) combining the first and second relationships to establish a system of equations 
in terms of the values at the distinct time points; and 

(f) finding responses of the circuit at the distinct time points by solving the 
1 5 established system of equations. 

2. The method of claim 1 , wherein each of the reference time points is one signal 
period away from its respective distinct point. 

20 3. The method of claim 1 , step (c) further comprising the steps of: 

selecting a set of initial values at the distinct time points; and 

performing integration from the distinct time points to their respective reference time 
points to generate values between the distinct time points and the reference time points. 

25 4. The method of Claim 3 further comprising a step of saving to disk Jacobian 

matrices at the integration time points. 

5. The method of claim 3, wherein the values generated by the integration are in 
time domain, step (c) further comprising the steps: 
30 generating values reflecting circuit responses in time domain based on the values 

generated by the integration. 
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6. The method of claim 1 , step (d) further comprising the steps of: 
selecting a set of initial values at the distinct time points; 

performing a forward multidimensional discrete Fourier transform on the initial values to 
transform the initial values from real values to complex values; 
5 performing a phrase shift of a sampling period on the complex values; and (12)- (1 3) 

performing an inverse multidimensional discrete Fourier transform to convert the shifted 
complex values to real values. 

7. The method of claim 1 , step (f) further comprising the step of: 

10 solving the established system of equations by performing a discrete numeric solution 

method to generate the circuit responses at the distinct time points. 

8. The method of claim 7 wherein in step (f) the discrete numeric solution method is 
Newton method. 

15 

9. The method of claim 7, wherein the established system of equations in step (e) is 
a system of nonlinear equations, the discrete numeric solution method contains solutions to a 
system of linear equations by applying matrix implicit iterative linear solvers to the linear 
equations, step (f) further comprising the steps of: 

20 (g) selecting a set of estimated values to reflect estimated circuit responses at the 

distinct time points; 

(h) establishing a system of linear equations at the estimated values; 

(i) preconditioning the system of linear equations to improve convergence of 
solution to the system of linear equations; and 

25 (j) solving the system of linear equations to generate correction values to adjust 

the estimated values. 

10. The method of claim 9, wherein step (i) comprises using a combination of a delay 
matrix and a linear model of the state transition function to construct the preconditioner. 
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1 1 . The method of claim 9, further comprising the steps of: 

(k) using the adjusted estimated values as newly estimated values to reflect 
estimated circuit response at the distinct time points; and 

(1) repeating steps (g>(k) until the estimated values are a satisfactory solution to 
5 the system of nonlinear equations established in step (e). 

12. The method of claim 1 1 , wherein die estimated values are in time domain, the 
method further comprising the step of: 

(m) converting the estimated values into frequency domain. 

10 

1 3 . The method of claim 1 , wherein the first and second sets of values are voltages at 
the distinct time points and reference time points. 

14. An apparatus for simulating responses of a circuit, the circuit receiving a periodic 
sample signal and at least one information signal, the apparatus comprising: 

(a) means for selecting a set of evenly spaced distinct time points; 

(b) means for defining a set of reference time points, wherein each of the 
reference time points is associated with one of the distinct time points; 

(c) means for establishing a first set of relationships between the values at the 
distinct time points and the values at the reference time points; 

(d) means for establishing a second set of relationships between the values at the 
distinct time points and the values at the reference points; 

(e) means for combining the first and second relationships to establish a system of 
equations in terms of the values at the distinct time points; and 

(f) means for finding responses of the circuit at the distinct time points by solving 
the established system of equations. 

15. The apparatus of claim 14, wherein each of the reference time points is one signal 
period away from its respective distinct point. 

30 
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16. The apparatus of claim 14, further comprising: 

means for selecting a set of initial values at the distinct time points; and 

means for performing integration from the distinct time points to their respective 

reference time points to generate values between the distinct time points and the reference time 

points. 

17. The apparatus of claim 16, wherein the values generated by the integration are in 
time domain, the apparatus further comprising: 

means for generating values reflecting circuit responses in time domain based on the 
values generated by the integration. 

1 8. The method of claim 1 4, further comprising: 

means for selecting a set of initial values at the distinct time points; 
means for performing a forward multidimensional discrete Fourier transform on the initial values 
to transform the initial values from real values to complex values; 

means for performing a phrase shift of a sampling period on the complex values; and 
means for performing an inverse multidimensional discrete Fourier transform to convert the 
shifted complex values to real values. 

19. The apparatus of claim 14, further comprising: 

means for solving the established system of equations by performing a discrete numeric 
solution method to generate the circuit responses at the distinct time points. 

20. The apparatus of claim 19 wherein the discrete numeric solution method is 
Newton method. 

2 1 . The apparatus of claim 1 9 wherein the established system of equations is a system 
of nonlinear equations, the discrete numeric solution method contains solutions to a system of 
linear equations by applying implicit matrix iterative linear solvers to the linear equations, the 
apparatus further comprising: 

(g) means for selecting a set of estimated values to reflect estimated circuit 
responses at the distinct time points; 

(h) means for establishing a system of linear equations at the estimated values; 
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(i) means for preconditioning the system of linear equations to improve 
convergence of solution to the system of linear equations; and 

0) means for solving the system of linear equations to generate correction values 
to adjust the estimated values. 

22. The apparatus of claim 2 1 , wherein the means (i) for preconditioning the system 
comprises means for combining the delay matrix and a linear model of the state transistion 
function to construct the preconditioner. 

23. The apparatus of claim 22, further comprising: 
(k) means for using the adjusted estimated values as newly estimated values to 

reflect estimated circuit response at the distinct time points; and 

(1) means for determining whether the estimated values are a satisfactory solution 
to the established system of nonlinear equations. 

24. The apparatus of claim 23, wherein the estimated values are in time domain, the 
method further comprising the step of: 

(m) means for converting the estimated values into frequency domain. 

20 25 - The apparatus of claim 14, wherein the first and second sets of values are voltages 

at the distinct time points and reference time points. 



15 



25 



22 



BNSOOCIO: <WO 0188830A1_U> 



WO 01/88830 



PCT/US01/15185 



1/8 



102 



is 
c 



104.1 





CIRCUIT 
WITH 14 NODES 


108 

/ 


104.2 

/ 




• 

• 
• 

104.S 

/ 


/UWWIAAAA/ 10 / 6 





FIG. 1 



BNSDOCID: <WO 0188830A1_I_> 



WO 01/88830 



PCT/US01/15185 




8NSDOCID: <WO. 



.018883QA1J_> 



WO 01/88830 



« 



PCT/US01/15185 




BNSDOCID: <W(D 



01B8830A1_L> 



WO 01/88830 



PCT/US0I/I5I85 



4/8 



2 
1.5 
1 

0.5 h 

0 
-0.5 

-1 
-1.5 



x BEFORE PRECONDITIONING 
o AFTER PRECONDITIONING 



-1.5 



10' 



= 10 



-1 



io 2 t: 



10 



r3 



j i i_ 



X 




-0.5 0 0.5 

FIG. 4 



i r- 



PRECONDITIONED 

RAW 



-i 1 1 1 1 1 i 1_ 



1.5 



0 50 100 150 200 250 300 350 400 450 500 550 

FIG. 5 



BNSOOCID: <WO 0188830A1 J_> 



WO 01/88830 



PCT/US01/15185 



CO 
-o 



20 
0 
-20 
-40 
-60 
-80 
-100 



-120 



-100 -50 



5/8 

i 1 1 r 



I 



-*— OUTPUT 
-e— INPUT 



50 100 150 200 250 300 350 400 
Hz 

FIG. 6 




59.97 59.98 59.99 



60 60.01 
MHz 

FIG. 7 



60.02 S0.03 60.04 



BNSDOCID: <WO O188830A1_l_> 



WO 01/88830 



6/8 



PCT/US01/15185 



( START 



802 



SELECT A SET OF DISTINCTIVE 
TIME POINTS 






DEFINING A SET OF REFERENCE 
TIME POINTS 






ESTABLISHING A FIRST SET OF 
RELATIONS 


i 




ESTABLISHING A SECOND SET OF 
RELATIONS 


1 




COMBINE THE FIRST AND SECOND SETS 
OF RELATIONSHIPS TO ESTABLISH 
A SYSTEM 






FIND RESPONSES ( 
SOLVING THE EST/ 


DF THE CIRCUIT BY 
\BLISHED SYSTEM 



804 



806 



808 



■810 



812 



814 



( END ^)-^" 



816 



FIG. 8 



BNSOOCID: <WO 01B883QA1_!_> 



WO 01/88830 



PCT/US01/15185 



7/8 



( START y ^ 902 



SELECT A SET OF ESTIMATED 
VALUES 



ESTABLISH A LINEAR SYSTEM 



PRECONDITIONING THE 
LINEAR SYSTEM 



SOLVE THE LINEAR SYSTEM 



ADJUST THE ESTIMATED VALUES 



916 

IS THE " 
SOLUTION SATISFACTORY^ 



904 



906 



908 



•910 



•912 



CONVERT THE ESTIMATED VALUES 
INTO FREQUENCY DOMAIN 



916 



( end 



918 



FIG. 9 



BNSDOCID: <WO 0188830A1_I_: 



WO 01/88830 



PCT/US01/15185 




BNSOOCID: <WO 0188830A1_I_> 



INTERNATIONAL SEARCH REPORT 



lot tional application No. 
PCT/US01/15185 



A. CLASSIFICATION OF SUBJECT MATTER 
IPC(7) :G06G 7/48; G06F 7/60 

US CL :703/2 f 3. 4, 5, 13. 14. 19, 20. 21. 22; 716/1. 2. 4 
According to International Patent Classification (IPC) or to both national classification and IPC 



B. FIELDS SEARCHED 



Minimum documentation searched (classification system followed by classification symbols) 
U.S.: 703/2. 3. 4. 5. 13, 14, 19. 20, 21 ,22; 71671. 2. 4 



Documentation searched other than niimmum documentation to the extent that such documents are included in the fields searched 



Electronic data base consulted during the international search (name of data base and. where practicable, search terms used) 



EAST. Proquest, IEEE. ACM 
search terms: mixed frequency-time, MPT, simulat* or model*, circuit* 



C DOCUMENTS CONSIDERED TO BE RELEVANT 



Category* 



Citation of document, with indication, where appropriate, of the relevant passages 



Relevant to claim No. 



A 
A 
A 



FENG et al. Efficient Computation of -Quasi-Periodic Circuit 
Operating Conditions via a Mixed Frequency/Time Approach. IEEE 
Proceedings of the 1999 Design Automation Conference. June 1999. 
pages 635-640; see entire document. 

US 5,995,733 A (ROYCHOWDHURY) 30 November 1999. 
US 5,588,142 A (SHARRIT) 24 December 1996. 

CHEN et al. A Mixed Frequency-Time Approach for Quasi-Periodic 
Steady-State Simulation of Multi-Level Modeled Circuits. 
Proceedings of the 1999 IEEE International Symposium on Circuits 
and Systems. May 1999. Vol. 6. pages 194-197. 



1-25 



1-25 
1-25 
1-25 



r j Father documents are listed in the continuation of Box C. Q See patent family annex. 



'O- 

•p* 



Special categories of cited documents: 

document defining the general state of the art which is not considered 
to be of particular relevance 

earlier document published on or after the international filing data 

document which may throw doubti on priority claim (s) or which is 
cited to establish the publication date of another citation or other 
special reason (as specified) 

document referring to an oral disclosure, use, exhibition or other 
means 

document published prior to the international filing date but later than 
the priority data claimed 



later document published after the international filing date or priority 
date and not in conflict with the application but cited to understand 
the principle or theory underlying the invention 

document of particular relevance; the claimed invention cannot be 
considered novel or cannot be considered to involve an inventive step 
when (he document is taken alone 

document of particular relevance; (he claimed invention cannot be 
considered to involve an inventive step when the document is 
combined with one or more other such documents, such combination 
being obvious to a person skilled in the art 

document member of the same patent family 



Date of the actual completion of the international search 



02 AUGUST 2001 



Name and mailing address of the ISA/US 
Commissioner of Patents and Trademarks 
Box PCT 

Washington, D.C 20231 
Facsimile No. -(703) 305-3230 



Date of mailing of< 



report 



Authorized officer 

SAMUEL BRODA 
Telephone No. <703) 305-1026 



Form PCT/IS A/210 (second sheet) (July 1998) * 



8MSOOCI0: <WO 0188830A1_I_> 



This Page is Inserted by IFW Indexing and Scanning 
Operations and is not part of the Official Record 



Defective images within this document are accurate representations of the original 
documents submitted by the applicant. 

Defects in the images include but are not limited to the items checked: 

□ BLACK BORDERS 

□ IMAGE CUT OFF AT TOP, BOTTOM OR SIDES 

□ FADED TEXT OR DRAWING 

□ BLURRED OR ILLEGIBLE TEXT OR DRAWING 

□ SKEWED/SLANTED IMAGES 



□ LINES OR MARKS ON ORIGINAL DOCUMENT 

□ REFERENCE(S) OR EXHIBIT(S) SUBMITTED ARE POOR QUALITY 

□ OTHER: 

IMAGES ARE BEST AVAILABLE COPY. 
As rescanning these documents will not correct the image 
problems checked, please do not report these problems to 
the IFW Image Problem Mailbox. 



BEST AVAILABLE IMAGES 




COLOR OR BLACK AND WHITE PHOTOGRAPHS 



Q.GRAY SCALE DOCUMENTS 



